%
% @brief 射频信号分析处理 相关的工具类
% @author 41204@qq.com / 2025-09-29
% @note 本类中的函数均以 Static 方式定义，以方便使用。
%

classdef rfsa    
    methods(Static)
        %% 绘制时域波形
        % @param iq_data [1xn complex double] IQ 数据
        % @param fs [number] 采样速率(Hz)
        function waveform(iq_data, fs)
            N = length(iq_data);  % 数据点数
            t = 1:N;
            t = t * 1e6 / fs;
            
            % 绘制时域波形
            figure;
            plot(t, 10*log10(abs(iq_data)));
            title('时域波形');
            xlabel('时间 (us)');
            ylabel('幅度 (dB)');
            grid on;
        end

        %% 频谱分析
        % @param iq_data [1xn complex double] IQ 数据
        % @param fs [number] 采样速率(Hz)
        % @param fshift [logical] 是否将频谱零频点移到中间位置, 默认为 false
        function [f, magnitude] = spectrum(iq_data, fs, fshift)
            if nargin < 3
                fshift = false; % shift 默认值
            end
            
            N = length(iq_data);  % 数据点数
            
            % 进行快速傅里叶变换（FFT）
            fft_data = fft(iq_data);
            clear iq_data

            if fshift
                % 调整频率轴
                f = fs * (-N/2:(N/2-1)) / N;
                mytitle = '频谱分析（频率轴居中）';
                % 将频谱零频点移到中间位置
                fft_data = fftshift(fft_data);
            else
                % 计算频率轴
                f = fs * (0:(N/2)) / N;  % 单边频率范围
                mytitle = '频谱分析（单边）';
                % 取正频率部分                
                fft_data = fft_data(1:N/2+1);  
            end

            % 计算幅度谱
            magnitude = abs(fft_data);            
            % 如果需要，计算功率谱
            %power = magnitude.^2 / N;
            
            % 绘制频谱
            figure;
            plot(f, 20*log10(magnitude));
            %plot(f, power);
            title(mytitle);
            xlabel('频率 (Hz)');
            ylabel('幅度 (dB)');
            grid on;
        end

        %% 计算射频接收机的幅度和相位响应
	    % 输入：
	    %   A: 输入IQ信号 [I_in, Q_in] 或 complex(I_in + 1j*Q_in)
	    %   B: 输出IQ信号 [I_out, Q_out] 或 complex(I_out + 1j*Q_out)
	    % 输出：
        %   h_freq: 复数形式的响应
		function h_freq = ap_response(A, B)
		    % 确保输入是复数形式
		    if size(A, 2) == 2  % 如果是分离的I,Q数据
			    A_complex = A(:,1) + 1j*A(:,2);
			    B_complex = B(:,1) + 1j*B(:,2);
		    else
			    A_complex = A;
			    B_complex = B;
		    end
    
		    % 计算FFT点数
		    N = length(A_complex);
    
		    % 计算输入输出的频谱
		    A_spectrum = fftshift(fft(A_complex, N));
		    B_spectrum = fftshift(fft(B_complex, N));
    
		    % 计算频率响应
		    h_freq = B_spectrum ./ A_spectrum;    
		    % 避免除零错误，将无穷大值替换为NaN
		    h_freq(isinf(h_freq)) = NaN;
		end
    
        % 绘制幅度和相位响应图
	    function plot_responses(h_freq, fs)    
		    figure;

            % 计算频率轴
            N = length(h_freq);
            freq = (-fs/2:fs/N:fs/2-fs/N)';

            % 计算幅度响应
            mag = abs(h_freq);
    
		    % 计算相位响应(弧度)
		    phase = angle(h_freq);
		    % 解卷绕相位
		    phase = unwrap(phase, 1.0*pi);
    
		    % 幅度响应
		    subplot(2,1,1);
            mag = 20*log10(abs(mag)); % 幅度响应转换为dB
		    plot(freq/1e6, mag, 'b', 'LineWidth', 1.5);
		    xlabel('频率 (MHz)');
		    ylabel('幅度响应 (dB)');
		    title('射频接收机幅度响应');
		    grid on;
		    xlim([-fs/2e6, fs/2e6]);
    
		    % 相位响应
		    subplot(2,1,2);
            phase = phase / (2*pi); % 相位响应转换为周数
		    plot(freq/1e6, phase, 'r', 'LineWidth', 1.5);
		    xlabel('频率 (MHz)');
		    ylabel('相位响应 (周)');
		    title('射频接收机相位响应');
		    grid on;
		    xlim([-fs/2e6, fs/2e6]);
		end
    end
end
